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1.  Introduction 


a.  Objective 

AFRL  has  the  mission  of  transitioning  a  next  generation  First  Principles  Predictive 
Assimilative  Operational  Satellite  Drag  Model  to  Air  Force  Space  Command.  This  Tech 
Memo  provides  an  initial  assessment  of  the  Thermosphere-Ionosphere  Electrodynamics 
General  Circulation  Model  (TIEGCM)  with  selected  supporting  scientific  research  noted 
as  needed  to  place  the  results  in  context;  it  is  not  intended  as  a  comprehensive  review 
paper.  Validation  is  done  with  satellite  accelerometer  neutral  density  data  from  the 
CHAMP  (-400  km),  GRACE  (-500  km)  and  historic  Air  Force  SETA  (-  200  km) 
missions.  Aspects  of  TIEGCM  Version  1.91  studied  are  solar  cycle,  day  of  year  and 
geomagnetic  storm  variations.  Also  described  are  assimilative  techniques  being 
implemented  for  the  latest  TIEGCM  model,  version  1.9.4.  We  identify  key  physics 
improvement  areas  for  next  generation  First  Principles  Forecast  operational  satellite  drag 
modeling  to  successfully  upgrade  current  and  future  satellite  orbital  and  reentry 
predictions. 

b.  Relevance 


Accurate  specification  and  prediction  of  the  space  environment  is  necessary  to  calculate 
satellite  drag  and  orbital  trajectories  precisely  to  support  Air  Force  Space  Command 
(AFSPC)  in  several  areas  including:  (1)  satellite  reentry  predictions,  (2)  timely  Collision 
Avoidance  warnings  needed  to  protect  International  Space  Station  and  other  space  assets 
from  being  hit  by  debris  and  other  space  objects  and  to  minimize  the  need  for  maneuvers 
and  (3)  catalog  maintenance  for  precision  orbit  determination  of  all  space  objects.  The 
AFSPC  present  stated  operational  goal  is  a  3 -day  satellite  density  forecast  with  error  of 
5%  for  altitudes  90  to  500  km.  As  a  more  quantitative  requirement,  Anderson,  et  al., 
2009  state  that  the  24  hour  orbit  prediction  errors  of  significance  include  250  meters  at 
200  km  and  100  meters  at  400  km.  Aerodynamic  drag  depends  on  neutral  density  p,  as 
well  as  area-to  mass  ratio  A/m,  drag  coefficient  CD  and  velocity  V  (includes  atmosphere 
rotation  with  earth  plus  residual  winds)  respectively: 


CDA 
V  m  ) 


pVA 


(1) 


Neutral  mass  density  variations  are  recognized  as  the  primary  cause  of  orbital  drag  errors. 
While  specific  knowledge  of  the  reentering  vehicle’s  area/mass  ratio  and  drag  coefficient 
(see  e.g.  Sutton,  2009;  Sentman,  1961),  as  well  as  neutral  wind  fields,  are  also 
critical,  we  focus  on  the  neutral  density  variations.  Thermospheric  density  is  driven 
mainly  by  the  highly  variable  sun’s  electromagnetic  radiation  and  interaction  of  solar 
particles  with  earth’s  magnetosphere.  Upward  propagating  waves  are  also  becoming 
identified  as  a  significant  contributor  to  many  aspects  of  thermosphere  variability.  Thus 
meaningful  progress  on  low  Earth  orbit  trajectory  forecasts  requires  significantly 
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improved  knowledge  of  space  weather  and  its  effects  of  Earth’s  neutral  upper 
atmosphere. 

c.  Thermospheric  Variability  Drivers 

The  Earth’s  upper  atmosphere  is  a  strongly  forced  and  coupled  system  driven  mainly  by 
the  sun’s  solar  EUV  and  solar  wind  interactions.  Solar  EUV  is  deposited  mainly  at  low 
and  middle  latitudes.  Geomagnetic  stonn  inputs  are  significant  at  high  latitudes.  Upward 
propagating  gravity  waves  and  tides  also  contribute  to  thermosphere  variability.  For  a 
review  of  fundamentals  of  thermosphere  heating,  see  Roble,  1977. 

Solar  EUV  is  the  main  mechanism  for  heating  Earth’s  thennosphere.  EUV  radiation,  at 
wavelengths  less  than  200  mn  typically  accounts  for  about  75  -  80%  of  the  energy  input 
to  the  thennosphere,  and  thereby  detennines  its  basic  structure.  The  neutral  density  varies 
by  an  order  of  magnitude  at  400  km  due  to  the  solar  variability  of  the  EUV  flux.  Below 
about  80  km  atmospheric  air  is  mainly  molecular  (~  78.08%  molecular  nitrogen,  20.95% 
molecular  oxygen  and  0.93%  argon,  with  trace  gases  including  0.038%  carbon  dioxide 
and  0.0005%  helium).  These  gases  are  mixed  such  that  the  composition  is  relatively 
constant  up  to  about  80  km.  Above  80  km  atomic  species  result  from  photodissociation 
and  chemical  reactions  due  to  solar  EUV  radiation.  At  higher  altitudes  the  different 
species  are  separated  gravitationally  such  that  the  main  thennospheric  constituent  is 
generally  atomic  oxygen  from  about  200  -  600  km.  Helium  densities  compete  with 
atomic  oxygen  above  about  600  km  during  low  solar  flux  conditions. 

Solar  EUV  fluxes  originate  in  the  chromosphere,  chromosphere-corona  transition  region 
and  corona.  In  contrast  to  visible  radiation  (i.e.  resulting  in  the  “solar  constant”),  EUV 
emissions  are  highly  variable,  with  chromospheric  emissions  changing  by  a  factor  of  2  or 
more  and  coronal  emissions  varying  by  a  factor  of  50-150  over  the  solar  cycle.  The 
thennosphere  response  to  solar  activity  consequently  varies  with  the  solar  cycle,  the  27 
day  solar  rotation  period,  with  day-to-day  changes  superimposed.  The  density  can  change 
by  an  order  of  magnitude  at  400  km  over  a  solar  cycle.  At  200  km  the  solar  cycle 
variability  is  ~  a  factor  of  3.  Solar  EUV  energy  is  deposited  mainly  at  low  to  mid¬ 
latitudes,  in  the  vicinity  of  the  sub  solar  point.  Since  this  heating  creates  a  pressure  bulge 
that  drives  winds  to  transport  heat  away  from  the  hot  dayside  toward  the  Earth’s  cold 
nightside,  temperatures  on  the  dayside  are  only  -30%  higher  than  those  on  the  nightside. 
Even  if  the  Earth  did  not  rotate,  solar  heating  would  still  reach  the  nightside  and  affect 
density  there. 

Geomagnetic  storms  account  on  average  for  about  20%  of  the  thennospheric 
heating.  These  storm  occurrences  are,  however,  episodic  and  unpredictable.  They  are  a 
result  of  solar  wind  energy  transferred  to  the  thennosphere  mainly  at  high  latitudes.  The 
solar  wind  is  a  supersonic  plasma  that  carries  with  it  the  magnetic  field  of  the  Sun. 
Interaction  of  the  solar  wind  with  Earth’s  magnetic  field  results  in  auroral  particle 
precipitation  and  frictional  dissipation  of  ionospheric  current  systems.  If  the  field 
direction  in  the  solar  wind  is  opposite  to  that  of  the  Earth’s  dipole  magnetic  field  then 
subsequent  interactions  allow  for  significant  amounts  of  mass,  momentum  and  energy  to 
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be  transferred  into  the  magnetosphere.  These  conditions  usually  result  in  magnetic  stonns 
that  produce  strong  auroral  particle  bombardment  and  strong  electric  currents  which  close 
in  the  ionosphere  and  drive  rapid  upper- atmospheric  motions.  During  geomagnetic 
storms,  the  disturbed  solar  wind  compresses  the  Earth’s  magnetosphere.  Intense  field- 
aligned  currents  couple  the  auroral  ionosphere  with  the  magnetosphere  extracting 
electromagnetic  energy  (Poynting  flux)  that  heats  both  the  ionized  and  neutral 
gases.  Joule  heating  heats  the  neutral  atmosphere  in  the  altitude  region  at  about  120-150 
km,  drives  horizontal  winds,  and  forces  upwelling  of  the  neutral  atmosphere.  Intense 
electric  fields  in  the  high-latitude  ionosphere  drive  rapid  plasma  convection  that  couple 
via  collisions  with  neutral  winds.  At  the  same  time,  the  auroral  oval  expands  and 
energetic  particles  precipitating  into  the  lower  thermosphere  enhance  ionospheric 
conductivities.  The  rate  of  electromagnetic  energy  input  at  high  latitudes  can  rise  to  more 
than  ten  times  greater  than  that  of  the  global  EUV.  The  resulting  heat  and  momentum 
input  produces  very  large  and  abrupt  changes  in  neutral  density  and  local  neutral 
composition,  generates  traveling  gravity  waves  and  excites  strong  winds.  The  heavier 
molecules  No  and  CL,  abundant  in  the  lower  thermosphere,  are  transported  to  higher 
altitudes.  Thus,  a  global  circulation  system  is  set  up  to  redistribute  mass,  momentum,  and 
energy  (Fuller-Rowell,  et  ah,  1996).  Localized  energy  input  is  spread  globally  with 
complex  temporal  and  spatial  variability. 

Figure  1  is  a  simplified  schematic  sketch  of  the  system  to  be  simulated  by  models, 
showing  the  relative  roles  of  high  latitude  solar  wind  heating;  low  latitude  solar  EUV  and 
upward  waves  in  various  latitude  regions  (Forbes,  2007).  Figure  2  (taken  from  Knipp,  et  ah, 
2005)  illustrates  the  relative  significance  of  solar  and  auroral  forcings  from  1975  to  2005. 
Solar  power  (grey  line)  generally  dominates  but  episodic  joule  power  during  geomagnetic 
storms  can  exceed  the  global  EUV  input.  Forcing  due  to  particle  precipitation  (grey  dots) 
is  relatively  small.  The  time  and  space  variation  of  neutral  atmosphere  density 
consequently  is  sensitive  to  solar  and  geomagnetic  activity  as  well  as  to  upward 
propagating  waves.  Most  thennospheric  variability  comes  from  variability  of  its  driving 
sources  over  periods  ranging  from  part  of  a  day  to  a  solar  cycle.  During  geomagnetic 
storms  especially,  the  thermosphere  can  vary  rapidly  and  irregularly  over  very  short  time 
periods. 

A  long-recognized  generally  less  significant  energy  source,  tidal  and  gravity  waves 
propagating  up  from  the  lower  atmosphere,  further  modulates  the  thermosphere.  Recent 
analyses  have  shown  that  upward  propagating  effects  originating  from  tropospheric 
heating  processes  that  reflect  land-sea  distributions,  topography,  and  other  factors  can 
potentially  explain  several  phenomena  including  longitudinal  density  variations  observed 
in  satellite  density  measurements  (e.g.  Forbes,  et  ah,  2009)  and  possibly  the  semiannual 
variation  (Akmaev,  2011). 

d.  Toward  Next  Generation  Satellite  Drag  Models 

Empirical  models  are  particularly  limited  in  their  capability  to  respond  to  the  highly 
variable  solar  and  geomagnetic  variations.  They  are  parameterized  in  terms  of  mainly 
proxy  geophysical  indices  that  are  daily,  3 -hourly,  and  hourly  at  their  highest  resolution, 
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and  so  are  not  suited  to  track  hour-to-hour  variations.  Further,  being  based  on  data,  they 
are  limited  when  extrapolating  beyond  the  range  of  data.  For  example,  the  Jacchia- 
Bowman  2008  (JB08)  operational  model  errors  for  large  storms  are  -21%  for  point 
measurements  (Bowman,  et  al.,  2008).  Using  orbit  averages  for  selected  stonns  showed  that 
the  smoothing  reduced  errors  to  -13%.  The  recent  solar  minimum  has  shown  errors  of  50% 
because  the  F10  -  EUV  one-to-one  relation  did  not  hold.  JB08,  while  using  measurements  as 
part  of  the  EUV  input,  is  still  heavily  weighted  in  favor  of  the  F10  index.  In  fact  the  model  has 
an  ad  hoc  correction  that  arbitrarily  represents  density  variability  with  decreasing  solar 
flux  in  order  to  better  match  data. 

The  scientific  and  operational  communities  have  recognized  that  significant 
improvements  will  require  a  shift  away  from  empirical  models  currently  used  to  physical 
models.  Thus,  next  generation  models  of  upper  atmospheric  density  forecasts  will 
necessarily  be  based  on  the  physics  of  solar  behavior  and  solar  terrestrial 
interactions.  Large  scale  computer  simulations  of  thermosphere  circulation  are 
demonstrating  their  ability  to  provide  global  distribution  of  the  mass  density  temperature 
and  winds.  These  complex  physical  models  hold  promise  of  great  progress  in  forecasting 
long  tenn  space  weather  variations.  While  the  climatological  variation  of  the  global  or 
orbit-averaged  density  is  relatively  well  captured  by  these  models,  describing  regional 
density  structures  and  short-tenn  density  variations  is  significantly  more  challenging.  The 
observational  solar  data  sets  and  underlying  subtle  physics  needed  to  achieve  the 
potential  accuracy  of  these  models  are  becoming  available.  The  satellite  drag  community 
recognizes  the  future  impact  of  these  models  vs  the  simpler  empirical  models  driven  by 
proxy  indicators  of  solar  heating. 

Physical  modeling  is  of  course  now  the  basis  for  our  improved  terrestrial  weather 
forecasts.  The  models  are  based  on  a  set  of  “first  principles”  equations  used  to  predict 
future  state  of  the  atmosphere:  equations  of  continuity,  conservation  of  energy, 
momentum  and  the  ideal  gas  law  are  used  to  evolve  the  thennospheric  parameters  neutral 
density,  temperature  winds  and  composition.  Their  use  is  obviously  not  straightforward 
in  the  empirical  modeling  sense  of  simple  analytical  functions.  The  nonlinear  partial 
differential  equations  are  necessarily  solved  by  numerical  methods  to  obtain  approximate 
solutions.  Weather  forecast  models  also  have  the  benefit  of  a  huge  amount  of  input  data 
from  land,  sea,  air  and  space-based  sensors  and  a  very  long  period  for 
upgrading.  Analogously,  space  weather  physical  models  use  the  same  basic  equations  to 
produce  space  environment  parameter  predictions  at  given  location.  As  noted  in  Section 
lc,  the  upper  atmosphere  is  strongly  affected  by  variable  high  energy  electromagnetic 
radiation  ultraviolet  from  the  Sun,  as  well  as  by  bombardment  from  energetic  auroral 
electrons  and  by  electric  current  flowing  from  the  magnetosphere.  It  is  also  influenced  by 
atmospheric  waves  that  propagate  up  from  lower  levels.  Unlike  the  lower  atmosphere,  the 
upper  atmosphere  is  heterogeneous,  with  diatomic  nitrogen  and  oxygen  giving  way  to 
monatomic  oxygen  at  the  higher  levels.  Also  unlike  the  lower  atmosphere,  the  dynamics 
is  strongly  influenced  by  molecular  viscosity,  heat  conduction,  and  diffusion,  as  well  as 
by  the  force  exerted  by  the  electric  current  flowing  through  the  Earth's  magnetic  field. 

Strong  external  forcing  and  relatively  rapid  dissipative  processes  mean  that  the  model 
cannot  run  freely  for  long  simulated  times  without  updates  to  the  inputs,  unlike  lower- 
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atmospheric  weather  and  climate  models.  Also  unlike  lower  atmosphere  models,  data- 
assimilation  techniques  are  challenged  by  sparse  thermosphere  and  heat  input  data. 

In  this  Tech  Memo  we  review  pertinent  aspects  of  neutral  density  validation  for  the 
NCAR  TIEGCM  (Thermosphere-Ionosphere  Electrodynamics  General  Circulation 
Model).  This  TM  uses  Version  1.91  of  the  TIEGCM.  Primary  references  for  the 
historical  development  of  this  model  are  Dickinson,  et  ah,  1981,  1984)  Roble  and  Ridley, 
1987,  1994,  Roble,  et  ah,  1987,  1988,  and  Richmond,  et  al.,  1992.  Another  excellent 
summary  tracing  the  model  development  to  its  current  version  1.9.4  is  given  by 
Solomon,  et  al.,  2012.  The  model  description  below  is  directly  from  NCAR’s  website  at 
http://www.hao.ucar.edu/modeling/tgcm/: 

“The  NCAR  Thermosphere-Ionosphere-  Electrodynamics  General  Circulation  Model 
(TIE-GCM)  is  a  comprehensive,  first-principles,  three-dimensional,  non-linear 
representation  of  the  coupled  thennosphere  and  ionosphere  system  that  includes  a  self- 
consistent  solution  of  the  low-latitude  electric  field.  The  model  solves  the  three- 
dimensional  momentum,  energy  and  continuity  equations  for  neutral  and  ion  species  at 
each  time  step,  using  a  semi-implicit,  fourth-order,  centered  finite  difference  scheme,  on 
each  pressure  surface  in  a  staggered  vertical  grid.  It  has  29  constant-pressure  levels  in  the 
vertical,  extending  from  approximately  97  km  to  500  km  in  intervals  of  one -half  scale 
height,  and  a  5°  x  5°  latitude-longitude  grid,  in  its  base  configuration.  The  time  step  is 
120  s.” 

“Hydrostatic  equilibrium,  constant  gravity,  steady-state  ion  and  electron  energy 
equations,  and  incompressibility  on  a  constant  pressure  surface,  are  assumed.  Ion 
velocities  are  derived  from  the  potential  field  created  by  combining  the  imposed 
magnetospheric  potential  with  the  low-latitude  dynamo  potential,  and  then  calculating  ion 
velocities  from  ExB  drifts,  rather  than  solving  the  ion  momentum  equations  explicitly. 
Some  minor  species  are  not  currently  included  in  the  model,  including  hydrogen  and 
helium  and  their  ions,  and  argon.  Several  parameterizations  are  used  in  the  TIE-GCM:  an 
empirical  model  is  used  to  specify  photoelectron  heating;  the  production  of  secondary 
electrons  is  included  using  an  empirical  model  derived  from  two-stream  calculations,  the 
effects  of  mixing  by  gravity  waves  are  included  using  an  eddy  diffusion  formulation; 

C02  is  included  by  specifying  a  lower  boundary  condition  and  assuming  that  it  is  in 
diffusive  equilibrium.  The  upper  boundary  conditions  for  electron  heat  transfer  and 
electron  number  flux  are  empirical  formulations.  At  the  lower  boundary,  atmospheric 
tides  are  specified  using  the  Global  Scale  Wave  Model  (GSWM).” 

AFRL  has  the  responsibility  of  Technology  Transitioning  to  AFSPC  an  operational 
physical  satellite  drag  model.  Figure  3  shows  a  block  diagram  of  an  assimilative 
TIEGCM.  In  practice,  an  operational  physical  model  uses  space-based  measurements 
of  actual  solar  EUV  and  high  latitude  heating  inputs  plus  tides  from  below  as  well  as 
thennosphere  and  ionosphere  data,  all  augmented  by  data  assimilation  techniques.  This 
effort  is  also  supported  by  the  five-year,  $1.5M/year  AFOSR  funded  MURI  with  the 
University  of  Colorado.  This  effort  completed  in  2012  significantly  advanced 
understanding  of  the  satellite  drag  environment  to  enable  specification  and  prediction  at 
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the  “next  level”  of  performance.  The  MURI  states  “  The  scope  of  this  MURI  topic  is  the 
elucidation  of  those  physical  concepts  and  predictable  key  indicators  of  energy  inputs  to 
the  atmosphere  that  will  update  and  calibrate  current  operational  satellite  drag  models  and 
lead  to  the  accurate  prediction  of  thermospheric  densities  and,  hence,  a  precise  prediction 
of  the  locations  of  satellites.”  Achieving  an  eventual  5%  forecast  capability  requires  (1) 
accurate  description  and  forecasts  of  relevant  heat  and  dynamics  inputs  (2)  a  model  that 
accurately  represents  the  thennosphere  response  to  the  external  heating  and  (3)  physical 
characteristics  of  the  spacecraft  object  used  in  orbit  propagators.  Figure  4  shows  how  the 
various  areas  of  study  under  the  MURI  support  the  capability  for  AFRL  to  transition  the 
needed  Precision  Orbit  Determination. 

e.  Memo  Overview 


This  memo  summarizes  the  initiation  of  detailed  evaluation  of  physical  satellite  drag 
models  to  support  implementation  of  next-generation  satellite  drag  forecast  capabilities. 
Section  2  briefly  describes  the  CHAMP,  GRACE  and  SETA  accelerometer  experiments 
that  provided  neutral  density  data  used  in  this  TM.  Sections  3,  4  and  5  illustrate  TIEGCM 
response  to  variations  with  day  of  year  (semiannual),  solar  cycle  and  geomagnetic  storms 
respectively.  Implementation  of  data  assimilation  techniques  are  noted  in  section  6.  A 
statistical  summary  of  results  and  preliminary  comparison  to  empirical  modeling  is  in 
Section  7.  A  summary  is  given  in  Section  8. 

2.  Data  Sources 

CHAMP,  launched  15  July  2000  into  a  near  circular,  near  polar  (inclination  angle  87°) 
orbit  with  an  initial  altitude  of  about  450  km,  provided  unprecedented  new  high- 
resolution  neutral  density  data  until  near  its  reentry  on  19  Sept  2011.  The  GRACE 
mission,  using  two  satellites  separated  by  a  nominal  220  km,  was  launched  17  March 
2002  into  an  89-degree  inclination  orbit  with  an  initial  perigee  near  500  km,  and  is  still 
operational.  CHAMP  uses  the  STAR  (Spatial  Triaxial  Accelerometer  for  Research) 
accelerometer  built  by  ONERA,  France  (Bruinsma,  et  al.,  2004).  GRACE,  at  higher 
altitude,  uses  SuperSTAR,  a  more  sensitive  version  of  STAR  (Tapley,  et  al.,  2004). 
Density  is  detennined  from  accelerometers  aboard  these  two  gravity  mission  satellites  by 
directly  measuring  decelerations  due  to  aerodynamic  drag.  Neutral  density  is  derived 
from  these  drag  measurements  according  to  Eq  (1)  using  knowledge  of  the  satellite’s 
area-to-mass  ratio,  drag  coefficient  and  velocity.  Data  processing  to  extract  neutral 
density  (and  winds)  from  these  accelerometers  is  given  in  Sutton,  et  al.,  2007. 

The  Satellite  Electrostatic  Triaxial  Accelerometer  instruments  were  flown  on  three 
satellites  in  near-circular,  sun-synchronous  orbit  (1030/2230  hours)  with  an  altitude  range 
between  170  and  240  km.  A  description  of  the  SETA  instrument  and  data  processing  to 
derive  density  and  winds  is  given  by  Marcos  and  Forbes,  1985  and  references  therein. 
Additional  results  from  the  SETA  flights  were  journal  published  between  1985  and  2000 
(see  e.g.  Rhoden,  et  al.,  2000  and  references  therein).  Figure  5  summarizes  the  satellite 
orbit  properties.  The  satellite  altitude  is  below  200  km  over  most  of  the  northern 
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hemisphere  dayside  conditions  throughout  each  mission.  We  utilize  data  at  latitudes  +/- 
83.7  degrees  from  the  SETA  2  (May-November,  1982). 

Data  are  presented  as  measured  and  model  densities  nonnalized  to  a  constant  altitude  of 
400  km  for  CHAMP  and  GRACE  and  to  200  km  for  SETA,  respectively.  In  both  cases 
nonnalization  is  carried  out  using  NRLMSIS  (appropriate  to  the  solar-geophysical 
conditions  for  the  measurement).  For  example  the  normalized  CHAMP  density  at  400 
km  is  calculated  as  CHAMP  (400  km)  =  CHAMP  (at  altitude  z)  times  the  ratio  of  MSIS 
(400km)/MSIS  (altitude  z).  The  effect  of  this  normalization  model  is  small,  particularly 
compared  to  the  spatial-temporal  variability  of  the  density  features  under  study.  In 
addition  ratio  of  measured  densities  to  TIEGCM  simulations  is  presented. 

3.  The  Semiannual  Variation 

The  most  prominent  feature  of  thennospheric  day  to  day  structure  (other  than 
solar/geomagnetic  variations)  is  the  so-called  semiannual  variation  (hereafter  SAV) 
characterized  by  maxima  near  equinoxes  and  minima  near  solstices.  It  was  first  observed 
in  satellite  drag  data  by  Paetzold  and  Zschomer,  1961.  Additional  studies  revealed  that 
the  magnitude  of  the  annual  minimum  to  maximum  variation  is  solar  cycle  dependent 
(Bowman,  et  ah,  2008)  and  can  be  more  than  100%.  The  maximum  yearly  difference,  from 
the  yearly  minimum  to  the  yearly  maximum,  varies  by  as  much  as  60%  from  year  to  year, 
and  the  phases  of  the  minima  and  maxima  also  can  change  by  about  a  month  from  year  to 
year.  There  are  also  latitude  structures,  hemispheric  asymmetries  (e.g.  Hedin  and  Mayr,  1987), 
and  interannual  variations  Over  50  years  since  the  SAV  discovery,  the  physical  mechanism  is 
not  clear.  A  full  discussion  is  beyond  the  scope  of  this  Memo  (see  Guo,  et  al.,  2008; 

Lei,  et  al.,  2012  for  recent  excellent  reviews). 

There  was  a  possibility  that  the  SAV  was  contained  in  existing  physical  models  (Roble, 
private  communication).  However,  Qian,  et  al.,  2009  utilizing  drag  data  from  two 
satellites  showed  that  the  SAV  was  not  contained  in  TIEGCM.  Figure  6  (taken  from 
Lin,  et  al.,  2010)  illustrates  this  model  deficiency.  The  top  part  of  the  figure  shows  ratio  of 
TIEGCM  to  density  data  derived  from  orbital  drag  on  two  different  satellites.  The  bottom 
part  of  the  figure  shows  ratio  of  TIEGCM  to  CHAMP  accelerometer  data  during 
2002.  In  both  cases  the  model  is  too  high  particularly  during  the  July  minima.  For 
CHAMP  the  difference  is  about  a  factor  of  two  higher  during  July  than  at  other 
times.  Note  that  the  data  also  showed  that  ratios  of  model  to  data  are  higher  during 
nighttime  than  daytime;  i.e.  there  is  an  apparent  local  time  error  in  the  model.  This  effect 
was  also  noted  in  the  assimilation  studies  described  in  Section  6.  Simulations  indicate 
that  seasonal  effects  in  TIEGCM  do  not  fully  account  for  the  observed  annual/semiannual 
amplitude,  primarily  because  of  the  lack  of  a  minimum  during  northern  hemisphere 
summer.  Qian,  et  al.,  2009  suggested  incorporating  a  variable  eddy  diffusion  coefficient 
larger  during  solstices  than  equinoxes,  and  stronger  turbulence  in  summer  than  in  winter 
tailored  to  match  the  orbital  drag  data  from  five  satellites.  TIEGCM  version  1.91  utilizes 
a  constant  value.  Lin,  et  al.,  2010,  using  CHAMP  data,  also  showed  how  a  variable  eddy 
diffusion  coefficient,  based  on  CHAMP  data  could  allow  TIEGCM  to  match  the 
semiannual  variation  in  CHAMP  data.  Their  variable  eddy  diffusion  coefficient  is  shown 
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to  the  left  on  Figure  7.  Increased  eddy  diffusivity  at  solstices  decrease  model  density  and 
decreased  eddy  diffusivity  at  equinoxes  will  increase  model  density.  As  a  result,  the  top 
right  panel  shows  CHAMP  (black  curve)  and  TIEGCM  (red  curve)  densities  during  2002 
now  in  good  agreement,  resulting  in  a  ratio  of  model  to  data  that  is  relatively  flat  (bottom 
right  panel).  All  TIEGCM  data  were  normalized  by  a  factor  of  1.15.  Seasonally  varying 
formulae  of  CHAMP  equatorial  eddy  diffusivity  with  a  systematic  bias  improves 
agreement  between  TIEGCM  and  CHAMP  neutral  density  values.  An  ad  hoc  prescribed 
seasonally  varying  eddy  diffusivity  correction  has  been  incorporated  into  version  1.9.4  of 
TIEGCM  (Qian  and  Solomon,  2012).  However,  as  was  noted  by  Qian,  et  ah,  2009,  this 
isn’t  necessarily  the  correct  or  complete  mechanism.  Upward  propagating  gravity  waves 
and  tides  may  be  involved.  The  seasonal  dependence  of  gravity  wave  activity  is  known  to 
peak  around  solstices,  and  minimize  around  equinoxes.  Akmaev,  et  ah,  2008  suggests 
gravity  wave  activity  in  the  Whole  Atmosphere  Model  (Akmaev,  2011)  may  account  for 
SAV. 

While  semiannual  variations  are  large,  the  changes  occur  over  periods  of  months.  Fitting 
the  periodicity  will  minimize  errors  in  forecasting  satellite  drag.  Hedin  and  Mayr,  1987 
in  an  analysis  of  wave  structures  in  Dynamics  Explorer  data  estimated  that  the  upward 
propagating  waves  contributed  about  5%  to  density  errors  at  low  latitudes.  However,  if 
the  mechanism  due  to  upward  propagating  effects  has  short  term  variations 
superimposed,  as  would  be  the  case  with  disturbances  originating  below,  then  these  will 
introduce  errors  into  the  density  forecasts.  Until  the  full  mechanism  is  understood  and  a 
system  exits  for  monitoring  it,  the  SAV  contribution  to  density  errors  remains  to  be 
resolved. 

4.  Solar  Cycle  Variations 

Annual  average  of  densities  normalized  to  400  km  densities  were  calculated  over  the 
period  2001-2009.  For  this  study  TIEGCM  was  driven  by  F 10  as  the  solar  EUV  input. 
Note  that  annual  solar  flux  varied  from  —181  in  2001  to  ~69  in  2008  and  2009.  To  focus 
on  model  perfonnance  without  contributions  due  to  large  geomagnetic  storms,  data  for  ap 
<_15  were  sorted  according  to  LT  06-18  (Blue)  and  18-06  hours  (Red).  Two  mean  values 
averaged  over  all  latitudes  were  determined  per  orbit  according  to  the  equatorial  LT.  The 
half  orbits  are  reasonably  consistent  with  algorithms  that  compute  orbit  averaged  data, 
and  were  also  used  in  a  preliminary  examination  of  day-night  differences.  Statistical 
analyses  over  the  solar  cycle  were  computed  based  on  the  half  orbit  averages. 

Plotted  in  Figure  8a  are  mean  value  ratios  of  TIEGCM/CHAMP  sorted  by  LT  bins  and  an 
average  NRLMSIS/CHAMP  (purple)  with  no  LT  bins.  The  TIEGCM/CHAMP  ratio 
shows  a  local  time  difference  with  night  ratios  higher  than  day  ratios  by  about  15% 
throughout  the  solar  cycle.  This  effect  was  noted  in  the  Section  3  above.  The 
TIEGCM/CHAMP  ratios  are  high  at  low  solar  flux,  as  expected,  because  densities  are 
anomalously  low  during  this  unusual  solar  minimum  (e.g.  Ernmert,  et  ah,  2010).  Since 
Figure  8a  shows  model  to  data,  ratios  are  higher  rather  than  lower.  The  EUV  input  is 
overestimated  by  the  F10  index  during  this  solar  minimum  so  model  values  do  not  catch 
this  anomaly.  Empirical  models,  also  driven  by  F10,  clearly  show  similar  behavior  (see 
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the  NRLMSIS  /Drag  curve).  While  part  of  this  decreased  neutral  density  could  be 
unmodeled  effect  of  thennospheric  cooling  related  to  anthropogenic  gases,  these  would 
be  of  the  order  20%  (Emmert,  et  ah,  2004;  Marcos,  et  ah,  2006).  Figure  8a  data  show 
differences  of  order  of  50%.  Therefore  the  F10  index  must  be  more  significantly 
underestimating  the  actual  EUV  flux  during  this  cycle,  i.e.  actual  EUV  flux  during  this 
cycle  is  lower  than  previous  ones  for  the  same  F10  values  (see  Solomon,  et  al.,  2010). 

Figure  8b  shows  relative  standard  deviations  (std  dev/mean)  over  the  solar 
cycle.  TIEGCM  errors  are  generally  higher  on  the  nightside  except  for  high  solar  flux 
(2002)  conditions.  TIEGCM  errors  are  typically  5-7%  higher  than  those  of  the  empirical 
NRLMSIS  model.  As  noted  at  the  end  of  Section  3  this  evaluation  is  based  on  version 
1.91  of  TIEGCM  and  therefore  the  error  is  increased  due  to  the  model’s  inaccurate 
semiannual  variation.  This  shortcoming  has  been  largely  overcome  in  newer  model 
versions;  it  is  reasonable  to  assume  that  the  TIEGCM  and  NRLMSIS  errors  are  more 
comparable.  An  updated  validation  is  recommended. 

The  limitations  of  using  F10  have  been  understood  since  the  earliest  empirical  models 
(Jacchia,  1964;  Jacchia,  1970).  Figure  8c  (left  side)  shows  solar  irradiance  vs  wavelength 
and  the  non-black  body  EUV  spectrum.  Note  that  the  F10  proxy  is  about  2  orders  of  magnitude 
off  the  right  hand  of  the  wavelength  scale.  The  bottom  half  of  this  chart  shows  ratio  of  EUV 
lines  as  solar  maximum  to  solar  minimum.  Some  EUV  lines  vary  by  an  order  of 
magnitude  over  the  solar  cycle  (Lean,  1997).  F10BAR  varies  by  about  a  factor  of  three 
over  the  solar  cycle.  Thus,  real  measurements  are  essential  for  capturing  thermosphere 
response  to  solar  electromagnetic  energy  input.  TIEGCM  can  accept  EUV  and  UV  data. 

Even  the  original  TGCM  was  developed  to  handle  57  EUV  lines.  While  routine 
parameterization  of  Solar  EUV  Experiment  (SEE)  data  as  a  physical  model  output  is  still 
in  a  research  state  (Dudock  deWit,  et  al.,  2009),  preliminary  studies  of  drag  data  from 
2003  have  demonstrated  the  benefit  of  realistic  SEE  inputs  in  TIEGCM  to  improve 
satellite  drag  specification  (Qian,  et  al.,  2008).  Further  analyses  of  drag  data  from  2002 
through  2007  by  Qian,  et  al.,  2009  found  that  TIEGCM  driven  with  measured  solar 
irradiances  does  an  excellent  job  in  reproducing  the  densities  forced  by  solar  cycle,  solar 
rotational,  and  geomagnetic  variation.  Solomon,  et  al.,  2010  used  SOHO  and  TIMED 
solar  data  to  show  that  EUV  levels  were  lower  in  2008  than  they  were  during  the 
previous  solar  minimum;  evidence  that  the  unusually  low  level  of  EUV  was  primary 
cause  of  the  observed  low  thermospheric  density.  Deng,  et  al.,  2012  noting  the 
anomalously  low  amount  of  geomagnetic  activity  in  2008  used  TIEGCM  simulations  to 
find  that  the  solar  irradiance  and  geomagnetic  energy  variations  account  for  3/4  and  1/4 
of  the  total  neutral  density  decrease,  respectively.  With  considerable  research,  more 
comprehensive  solar  EUV  input-drag  response  relations  will  become  achievable  using 
new  observational  tools  including  the  EVE  experiment  (Woods,  et  al.,  2012)  on  Solar 
Dynamics  Observatory  plus  abundant  satellite  drag  and  satellite  accelerometer  density 
data. 
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5.  Geomagnetic  Storm  Response 


5a.  Section  overview 


Thermosphere  response  to  geomagnetic  conditions  involves  a  series  of  complex  and 
interrelated  processes  with  high  latitude  solar  wind  deposition  of  energy  being  converted 
to  global  neutral  density  perturbations  as  noted  in  Section  lb.  The  serious  shortcomings 
of  present  empirical  models  are  most  evident  during  large  magnetic  storms  when  these 
complex  inputs  and  thennospheric  responses  become  difficult  to  capture.  Geomagnetic 
activity  is  recognized  as  the  largest  satellite  drag  uncertainty:  AFSPC  has  shown  that 
tracking  errors  at  400  km  during  even  moderately  disturbed  periods  are  65%  greater  than 
at  quiet  times.  In  this  section  we  illustrate  general  features  of  TIEGCM’s  ability  to 
represent  the  neutral  density  response  to  geomagnetic  activity  and  use  statistical 
evaluations  to  point  out  current  capabilities  and  areas  of  improvement  under 
investigation.  We  use  the  higher  altitude  CHAMP  and  GRACE  data  to  demonstrate  storm 
responses  at  high  and  low  solar  activity,  then  show  SETA  1982  data  to  estimate  TIEGCM 
capabilities  in  the  reentry  region.  Figure  9  shows  CHAMP  (blue  line)  &  GRACE  (red 
line)  perigee  altitude  vs  time  for  2001  through  2008  (top  chart)  with  corresponding  daily 
F10  solar  flux  values  (bottom  chart).  Selected  storm  periods,  noted  by  the  vertical  blue 
lines  are  in  2001-2  (F10-150  to  >200  sfu)  and  2007-2008  (F10  ~65  sfu).  All  storm 
periods  are  near  November  to  reduce  seasonal  differences  in  the  solar  flux  comparisons. 
For  high  solar  flux  there  are  two  cases  for  CHAMP  in  2001  with  ap  of  48  and  39 
respectively  and  one  in  2002  for  CHAMP  and  GRACE  with  ap  =  48.  For  solar  minimum, 
there  is  one  CHAMP-GRACE  case  in  2007  with  ap  =  56  and  one  in  Nov  2008  with  ap 
=32. 

5b.  Storm  Effects  at  High  Altitudes  During  High  Solar  Flux 

Figure  10  shows  in  geographic  latitude-time  coordinates:  (a)  CHAMP  density  response 
for  two  moderate  storms  in  2001  (FBAR  -214),  (b)  TIEGCM  simulations  of  these  storms 
and  (c)  ratio  of  CHAMP  data  to  TIEGCM  model.  For  the  stonn  period  starting  Day  304; 
the  top  left  side  of  Figure  10a  shows  the  early  evening  hemisphere  of  the  orbit  (1924 
hours)  and  the  chart  below  it  (Middle)  gives  the  early  morning  side  of  the  orbit  (0724 
hours).  The  normalized  densities  are  obtained  by  first  normalizing  all  CHAMP  densities 
to  400  km  altitude,  then  comparing  stonn  period  densities  to  a  “quiet”  condition 
detennined  from  pre-storm  density  values  vs  latitude  for  the  previous  day  where  ap  was 
low  and  fairly  constant.  The  bottom  charts  show  ap  values  for  each  case  are 
similar:  maximum  ap=  48  on  day  305  (left)  and  ap  =  39  on  day  358  (right).  In  both 
cases,  the  maximum  density  increase  is  at  high  latitudes  (up  to  factor  of  two)  for  both 
local  times  and  there  is  greater  latitude  penetration  of  storm  enhancements  on  nightside. 
The  response  to  the  second  stonn  starting  on  day  357  is  much  weaker  (right  side  of 
Figure  10a),  indicating  possibly  a  local  time  dependence  as  well  as  the  relatively 
imprecise  correlation  of  heat  input  simulation  vs  actual  drivers.  The  TIEGCCM 
simulated  response,  Figure  10b,  shows  similar  features  but  the  response  is  weaker.  Ratios 
of  CHAMP  2001  densities  to  TIEGCM  are  shown  in  Figure  10c.  The  large  amounts  of 
blue  on  the  chart  indicate  that  the  model  has  the  relative  responses  fairly  well  represented 
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particularly  early  afternoon  but  with  some  underestimation  for  both  storms  at  high 
latitude  -dawn/  and  dusk  where  there  are  localized  errors  of  a  factor  of  two.  The 
equatorial  penetration  is  also  under  represented  at  night  and  early  morning  hours. 

We  next  compare  the  response  to  a  geomagnetic  storm  for  both  CHAMP  and  GRACE  in 
Figure  1  la.  During  the  period  studied,  days  313-316  in  2002,  CHAMP  is  in  a  0917/21 17 
hours  orbit  (top  left  and  right  respectively)  and  GRACE  is  in  an  0815/0615  orbit  (bottom 
left  and  right  respectively).  Normalized  density  responses  (as  for  Figure  10a)  for  this 
moderate  storm  during  high  solar  flux  conditions  (FBAR  =  167;  ap  max  =  48  on  day  314) 
show  similar  features  for  CHAMP  and  GRACE.  As  for  the  CHAMP  2001  data  in  Figure 
10,  the  maximum  density  increase  is  at  high  latitudes  (up  to  factor  of  two)  for  all  local 
times.  Greater  latitude  penetration  of  storm  enhancements  occurs  on  nightside.  Figure 
lib,  ratios  of  both  CHAMP  and  GRACE  to  TIEGCM,  show  model  simulations 
underestimate  response  to  this  stonn,  with  the  largest  errors,  up  to  a  factor  of  two  as  for 
CHAMP  in  2001,  mainly  confined  to  the  high  latitude  winter  hemisphere. 

5c.  Storm  Effects  at  High  Altitudes  During  Eow  Solar  Flux 

The  response  to  a  moderate  storm  in  2007  (FI 0=68,  FBAR  =  72)  with  max  ap  =  56  on 
day  324  (20  Nov)  is  shown  in  Figure  12a  for  CHAMP  (top  left  at  0826  hours  and  top 
right  at  2026  hours)  and  GRACE  below  at  0123  hours  (left)  and  1323  hours  (right).  There 
are  altitude  and  local  time  differences  in  the  responses  seen  by  the  two  satellites.  As  a 
general  feature,  at  these  low  solar  flux  conditions,  large  density  increases  are  up  to  factor 
of  three  (compared  to  -2  at  high  solar  flux)  during  the  approximately  day  long 
perturbation.  The  response  at  GRACE  altitudes  is  overall  stronger  than  that  for  CHAMP. 
This  may  be  partly  attributed  to  different  composition  effects.  Thayer,  et  ah,  2012  note 
the  impact  of  Helium  on  geomagnetic  response  at  GRACE  altitude.  TIEGCM  does  not 
include  Helium.  Corresponding  simulated  TIEGCM  densities  in  Figure  12b  show 
enhanced  response  in  the  northern  (summer)  hemisphere  extending  to  all  latitudes,  but  a 
weaker  extension  for  CHAMP  at  0826  hrs.  Ratios  of  CHAMP  and  GRACE  2007 
densities  to  TIEGCM  (Figure  12c)  show  that  the  model  again  tends  to  underestimate 
responses  at  all  latitudes  and  at  all  local  times  by  as  much  as  a  factor  of  two. 

The  solar  cycle  response  study  above  revealed  a  greater  relative  response  magnitude  at 
solar  minimum  attributed  to  similar  energy  absorbed  by  less  mass  (Burns,  et  ah,  2004)  and 
lower  scale  height  and  to  increased  disturbance  with  altitude  due  to  lower  scale  heights 
during  solar  minimum.  The  disturbance  region  extension  to  all  latitudes  during  solar  min 
appears  related  to  greater  horizontal  advection.  Weaker  EUV-driven  pressure  gradients 
result  in  stronger  equatorward  winds  from  Joule  heating.  Also,  very  low  solar  min  density 
data  imply  weak  summer-winter  circulation,  also  enhancing  summer  hemisphere 
disturbance  equatorward  propagation.  It  is  noted  that  the  mean  density  is  lower  by  around 
an  order  of  magnitude  so  the  satellite  drag  impact  is  not  a  significant  operational 
problem. 

Figure  13a  shows  CHAMP  &  GRACE  normalized  response  to  a  small  geomagnetic 
disturbance  (max  ap=32)  during  days  312-315  of  the  deep  2008  solar  minimum  (F10=69; 


Approved  for  public  release;  distribution  is  unlimited. 


11 


FBAR  =68). The  satellites  are  in  very  similar  local  time  orbits  near  roughly  1000/2300 
hours.  The  responses  are  not  readily  interpreted:  the  maximum  responses  occur  on  the 
nightside  for  both  satellites  and  are  at  low  to  middle  latitudes.  Further  study  is  needed  to 
unravel  this  low-mid  latitude  increase  observed  by  both  satellites.  On  the  dayside,  the 
response  magnitude  tends  to  favor  northern  latitudes.  GRACE  response  is  higher  than 
that  of  CHAMP  at  lower  latitudes.  The  relative  response  is  again  by  a  factor  of  three.  The 
ratios  of  data  to  model  (Figure  13b)  are  again  up  to  a  factor  of  two. 

5d.  Storm  Effects  in  Reentry  Region 

We  next  examine  TIEGCM  storm  response  in  the  reentry  region.  SETA  data  remain  the 
only  source  of  continuous  full  orbit  accelerometer  data  in  this  region.  As  a  case  study, 

Figure  14a  shows  SETA-2  measured  density  nonnalized  to  an  altitude  of  200  km  vs 
latitude  and  time  during  the  August  1982  period  with  maximum  ap  of  179.  The  three 
frames  are:  dayside  density  (top),  nightside  density  (middle)  and  geomagnetic  indices  ap 
and  Dst  (bottom).  The  dayside  maximum  response  of  -65%  occurs  near  70  deg  north 
latitude.  In  the  south  (winter)  hemisphere  the  maximum  response  is  -30%,  presumed  due 
to  reduced  conductivity  and  hence  reduced  joule  heating.  The  response  penetrates  to 
lower  latitudes  in  the  summer  hemisphere  due  to  prevailing  summer  to  winter  hemisphere 
winds.  On  the  nightside  penetration  to  low  latitudes  is  about  the  same  in  both 
hemispheres  due  to  auroral  winds  reinforced  by  day-night  flow  and  reduced  ion  drag. 

Figure  14b,  Top,  gives  the  SETA  response  now  normalized  to  pre-storm  conditions  and 
again  showing  maximum  increase  at  high  latitudes  >50%.  and  greater  latitude 
penetration  of  storm  enhancements  on  nightside.  Ratio  of  measured  SETA  data  to 
TIEGCM  (both  at  altitude)  is  on  the  bottom  half  of  Figure  14c.  Here,  the  ISEE-3  E-field 
data  were  available  to  drive  the  Weimer  2005  Model  (REF).  Dayside  underestimates  are 
generally  -10%;  some  high  latitudes  errors  up  to  -40%.  Nightside  typical  errors  are  from 
-30%  up  to  50%. 

In  Figure  15  we  further  examine  this  storm  showing  density  and  model  data  at  three 
geographic  latitudes  (20,  60  and  70  deg)  in  the  dayside  northern  (left)  and  southern 
(right)  hemispheres.  Empirical  models  densities  are  included  for  comparison.  Measured 
density  at  200  km  (Y  axis)  has  more  structure  than  the  JB08,  J70,  NRLMSIS  and 
TIEGCM.  All  models  underestimate  the  response  at  higher  latitudes.  For  example,  at  60 
deg  in  the  northern  hemisphere,  the  ratio  of  density  for  day  219.5,  approximately  the 
maximum  time  of  maximum  response  to  this  storm,  to  day  218,  a  quiet  period,  is  36%  for 
SETA,  25%  for  NRLMSIS,  9%  and  16%  for  TIEGCM. 

6.  Assimilation  Techniques 

As  was  shown  in  Figures  3  and  4  data  assimilation  is  a  critical  part  of  operational  model 
development.  The  concept  of  dramatically  reducing  empirical  model  satellite  orbit 
prediction  errors  by  assimilating  drag  data,  developed  by  AFRL  (Marcos,  et  al.,  1998),  has 
been  very  successfully  demonstrated  in  AFSPC’s  HASDM  model  (Storz,  et  al.,  2002)  and  led 
to  the  MURI  program.  Prior  to  this,  physical  models  of  the  thermosphere  almost  exclusively  dealt 
with  neutral  composition  and  winds.  Now  they  deal  most  heavily  with  satellite  drag  data.  After 
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a  transition  period  of  several  years,  physical  models  with  assimilation  of  large  amounts  of 
observations,  are  responsible  for  our  local  weather  forecasting.  Our  next  generation 
Space  Weather  forecasting  is  dependent  on  developing  data  assimilation  techniques 
(Matsuo,  et  ah,  2012  and  references  therein).  Assimilative  modeling  provides  a  means  to 
systematically  identify  and  correct  the  inconsistencies  between  model  specification  and 
observations.  For  example,  as  shown  in  Section  5  above,  today’s  models  are  capable  of 
reproducing  generic  geomagnetic  storm  effects;  modeling  specific  storms  is  challenging 
due  to  lack  of  accurate  descriptions  of  the  storm  energy  inputs.  Assimilation  of 
measurements  into  the  modeling  permits  a  means  for  compensating  for  the  uncertainty  in 
model  inputs  for  a  given  period  and  addressing  the  most  pressing  satellite  drag  need — 
accurate  speciation  and  forecasts  of  density  variability  due  to  geomagnetic  activity. 
Kalman  Filter  (Kalman,  1960)  techniques  have  been  developed  to  optimally  combine  the 
understanding  captured  in  physical  models  with  available  data  to  specify  and  forecast  the 
state  of  the  system  (Matsuo,  et  ah,  2012  and  references  therein).  Current  thermospheric 
efforts  necessarily  involve  assimilating  limited  high  resolution  satellite  drag 
measurements  into  a  first  principles  model  to  achieve  more  accurate  global  neutral 
density  specification. 

Supported  by  AFRL’s  Orbital  Drag  Environment  program  an  assimilative  model  is  being 
developed  (Retterer,  2012)  for  an  operational  TIEGCM  with  a  unique  approach  using  an 
ensemble  Kalman  filter  to  ingest  satellite  drag  and  other  data.  The  goal  is  to  combine 
observational  data  and  theoretical  knowledge  to  obtain  the  more  complete  and  accurate 
thennosphere-ionosphere  description  needed  to  meet  Air  Force  requirements.  The  effort 
comprises  several  parts:  the  first -principles  thermospheric  and  ionospheric  models,  the 
assimilation  framework  and  the  ingested  data.  Details  will  be  provided  in  a  separate 
report. 

An  initial  test  was  done  using  the  most  recent  version  (1.9.4)  of  the  TIEGCM  model  over 
28  days  starting  on  day  290,  2004  (and  thus  including  the  big  storms  of  November 
2004).  Neutral  densities  were  compared  with  CHAMP  at  the  satellite  location.  The 
objective  was  examination  of  parameterizations  of  TIEGCM  inputs  for  the  assimilation 
process.  The  parameters  first  tried  include  1)  a  scale  factor  for  the  solar  EUV  flux 
(sflxfac),  2)  a  scale  factor  between  the  solar-wind  electric  field  and  the  Joule  heating 
(joulefac),  and  3)  a  factor  for  the  eddy  diffusivity  (eddyfac).  The  choice  of  the  first  factor 
sflxfac  is  obvious.  The  second  scale  is  chosen  because  the  Joule  heating  is  a  nonlinear 
function  of  the  electric  field,  which  can  be  readily  specified  using  higher  temporal 
resolution  (5  minute)  solar- wind  specification  files  from  the  ACE  data.  The  third  scale 
factor  is  selected  because  it  is  empirically  specified,  anyway.  The  idea  is  to  adjust  scale 
factors  for  various  physical  processes  to  capture  the  natural  variability  that  isn't 
represented  by  the  variation  of  the  parameters  used  to  drive  the  model.  Trials  were  made 
with  all  possible  combinations  of  4  values  of  sflxfac,  3  values  of  joulefac,  and  three 
values  of  eddyfac.  The  parameter  set  producing  the  minimum  error  was  selected  as  the 
optimal  parameter  choice  for  that  day.  This  approach  addresses  how  uncertain  model 
parameters  can  be  optimized  to,  on  average,  give  the  best  agreement  with  observations.  It 
will  allow  to  best  compare  model  variability  including  information  about  scale  sizes  with 
variability  in  observations  to  determine  model  fidelity.  Finally  it  will  improve  observed 


Approved  for  public  release;  distribution  is  unlimited. 


13 


spatial/temporal  model  input  parameters  to  get  complete  specification  of  the  boundary 
conditions. 

Initial  tests  using  the  parameters  above  found  that  this  assimilation  process  could  produce 
good  agreement  in  the  daytime  or  nighttime,  but  not  in  both  portions  of  the  local  day 
simultaneously.  If  the  daytime  densities  were  right,  then  the  nighttime  densities  would  be 
too  high.  Some  additional  parameter  was  required.  The  relevant  process  must  be  one 
that  reduces  the  density  when  there  is  no  solar  heating;  given  the  uncertainties  in  rates, 
we  adopted  the  NO  cooling  rate  as  an  additional  factor.  To  explore  the  effects  of 
adjusting  this  rate,  we  employed  the  measured  SEE  solar  spectrum  to  drive  TIEGCM 
instead  of  using  the  EUVAC  spectrum  parameterized  by  FI  0.7,  to  eliminate  one  of  the 
other  adjustable  factors. 

Results  for  the  accuracy  of  the  TIEGCM  model  in  reproducing  CHAMP  neutral  density 
data  when  the  density  data  were  used  to  set  simple  scaling  of  parameters  of  the  model  are 
summarized  in  Figure  16  a-d.  Figure  16a  shows  the  density  over  several  orbits  during  a 
portion  of  day  294,  2004,  a  geomagnetically  quiet  day  for  the  CHAMP  measurements 
and  several  unconstrained  models.  The  top  panel  gives  the  mass  density,  the  middle 
panel  the  satellite  latitude,  and  the  bottom  panel  the  longitude  (black  curve)  and  local 
time  (blue  curve)  of  the  satellite.  In  the  top  panel,  the  red  curve  is  the  mass  density 
observed  by  CHAMP,  the  black  curve  is  the  model  density  from  the  empirical  NRLMSIS 
model,  the  pink  curve  is  TIEGCM  using  the  measured  SEE  solar  spectrum,  while  the  blue 
and  green  curves  are  the  TIEGCM  using  the  EUVAC  solar  spectrum  parameterized  by 
FI 0.7.  The  model  runs  for  green  curve  used  the  standard  OMNI  database  for  geophysical 
parameters  (FI 0.7,  Kp,  ACE  solar- wind  parameters),  while  the  model  runs  for  the  blue 
curve  used  a  geophysical  database  at  a  finer  time  resolution.  This  figure  shows  that  the 
parameterized  solar  spectrum  EUVAC  drives  TIEGCM  with  noticeably  poorer  results 
than  the  measured  solar  spectrum.  TIEGCM  using  SEE  spectrum  produces  reasonable 
daytime  densities  but  a  reduced  day/night  dynamical  range. 

Figure  16b  shows  the  same  quantities  as  Figure  16a,  but  for  the  TIEGCM  runs 
constrained  by  adjusting  rate  factors  to  better  match  the  density  data.  The  black  and  pink 
curves  are  again  the  MSIS  and  unconstrained  TIEGCM/SEE  runs.  The  blue  curve  is  now 
the  TIEGCM/SEE  run  with  the  NO  cooling  factor  adjusted  to  improve  the  match  with 
densities.  This  fitting  required  an  increase  in  the  NO  cooling  rate  by  about  a  factor  of 
three;  each  day  was  adjusted  individually,  but  the  resulting  NO  cooling  rate  factor  was 
about  the  same  every  day.  For  the  run  shown  by  the  green  curve,  an  additional  factor,  the 
Joule  heating  rate  factor,  was  also  permitted  to  vary.  The  improvement  of  the  model  fit 
was,  however,  minimal. 

Figure  16c  shows  the  statistical  results  over  the  whole  30-day  trial  for  the  TIEGCM 
model  using  a  measured  solar  spectrum  (the  SEE  data),  so  there  was  no  degree  of 
freedom  for  adjusting  the  solar  radiation  intensity.  The  blue  points  show  the  results  when 
the  scaling  of  the  NO  cooling  rate  was  adjusted  to  give  the  best  fit;  the  green  points  show 
the  results  fitting  both  the  NO  cooling  rate  and  the  scale  factor  of  the  Joule  heating  rate, 
while  for  comparison  the  pink  points  show  the  results  for  the  TIEGGM  densities  when  no 
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adjustment  of  any  parameter  was  made,  and  the  black  points  show  the  results  for  an 
empirical  climatological  model,  MSIS.  The  top  panel  gives  the  daily  mean  density  along 
the  CHAMP  orbit  (red  points),  along  with  the  TIEGCM  model  evaluated  at  those  points. 
The  second  panel,  labeled  prediction  efficiency,  gives  the  difference  between  the  mean- 
square  variations  of  the  density  and  the  mean  square  error  of  the  model  calculation, 
nonnalized  by  the  mean  square  variation  of  the  density.  Thus,  if  this  quantity  were  unity, 
the  model  was  able  to  track  the  variations  of  the  measurements;  smaller  values  are  worse, 
and  if  the  quantity  is  negative,  then  the  errors  in  the  model  are  worse  than  the  natural 
fluctuations  in  the  observed  quantity.  This  plot  shows  that  MSIS  actually  performs  the 
best  most  of  the  time  in  predicting  the  variations  of  the  density,  but  fails  during  the  time 
of  geomagnetic  activity  (days  312-316  of  2004).  The  fitted  TIEGCM  model  perfonns 
adequately  during  the  quiet  days,  but  is  the  best  during  the  period  of  geomagnetic 
activity.  Fitting  the  Joule  heating  rate  improves  the  perfonnance  only  slightly  (the  green 
points  are  a  little  better  than  the  blue  ones);  the  2-parameter  fit  may  not  be  better  because 
of  the  fine  time  step  used,  implying  that  the  mean  values  describe  the  electric  fields  well. 
The  final  set  of  points  (pink),  are  for  the  TIEGCM  model  without  any  adjustment,  and 
they  are  the  worst  in  accuracy.  The  other  two  panels  of  the  plot  show  the  mean  and  root- 
mean-square  errors  of  the  different  model  runs. 

Finally,  Figure  16d  shows  a  sample  of  the  difference  in  neutral  winds  at  240  km  altitude 
between  two  model  runs  (fitted  and  unfitted)  during  an  interval  of  geomagnetic  activity. 
The  constrained  model  had  the  NO  cooling  factor  fitted  to  have  its  densities  better  match 
the  CHAMP  densities.  Differences  in  winds  over  much  of  the  globe  have  magnitudes  up 
to  200  m/s.  Because  satellite  drag  depends  on  the  velocity  of  the  satellite  relative  to  the 
ambient,  these  errors  in  winds  contribute  to  the  satellite  errors  we  want  to  reduce. 

7.  Statistical  Comparison  of  TIEGCM  and  Empirical  Models 

A  statistical  analysis  of  data  vs  TIEGCM  version  1.91  provides  a  more  quantitative 
estimate  of  the  model  capability  to  describe  thermosphere  variability.  We  use  CHAMP 
2002-2005  and  SETA-2  1982  point  density  values  and  focus  on  the  northern  hemisphere 
since  that  is  where  SETA  is  most  accurate.  Data  are  first  sorted  into  day  and  night  local 
time  bins;  6-18  hours  day  and  not  within  6-18  hours  for  night  for  CHAMP.  SETA  day 
and  night  bins  are  approximately  1030  and  2230  respectively.  Following  the  day  /night 
binning,  all  data  are  sorted  into  six  ap  bins:  0-10,  10-25,  25-50,  50-100,  100-200  and  200- 
400.  We  then  examine  the  low  ap  bin  and  the  two  highest  ap  bins  combined  (ap>100)  as  a 
function  of  latitude.  In  this  part  of  the  study  the  data  are  put  into  10  degree  geographic 
latitude  bins.  Then  a  mean  and  standard  deviation  are  calculated  for  each  bin.  We  then 
compared  CHAMP  and  SETA  data  to  TIEGCM,  JB08  and  other  models  that  were 
described  in  other  studies  (Marcos,  et  al.,  2009).  CHAMP  data  come  from  an  AFRL 
report  to  Rand  Corp  (Marcos,  et  al.,  2009);  SETA  data  are  from  Marcos  and  Forbes, 

1985.  This  provides  a  measure  of  how  well  the  physical  model  used  in  this  TM 
compares  to  the  current  operational  model.  Assimilation  results  were  not  used  in  this 
comparison.  The  comparison  may  generally  not  match  exactly  for  the  databases,  but 
differences  should  be  small. 
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Figure  17a  shows  CHAMP  to  TIEGCM  mean  ratios  (top)  and  standard  deviations 
(bottom)  for  day  (left)  and  night  (right)  local  time  bins.  Day  data  are  in  blue  and  night 
data  are  in  red.  Note  that  the  data  for  each  ap  bin  are  plotted  at  the  beginning  interval  of 
the  bin,  rather  than  in  the  middle.  The  mean  values  increase  with  increased  ap,  from 
about  1  in  the  lowest  bin  to  about  1.3  in  the  highest  bins,  indicating  the  observed 
underestimation  of  the  density  response  with  increasing  geomagnetic  activity.  Similarly, 
the  standard  deviations  increase  with  geomagnetic  activity  as  would  be  expected  from 
about  25%  for  the  0-10  ap  bin  to  about  3 1%  for  ap  >  100.  Dayside  means  tend  to  be  about 
10%  higher  than  night  means  apparently  reflecting  the  higher  nighttime  model  densities 
implied  in  semiannual  variation,  solar  cycle  and  assimilation  studies  described 
previously.  Standard  deviations  tend  to  be  about  4%  lower  during  daytime.  In  Figure  17b, 
for  comparison,  CHAMP  mean  ratios  to  JB08  are  similar;  being  about  1.07  in  the  0-10  ap 
bin  and  very  close  to  unity  above  ap  of  100.  Standard  deviations  increase  from  16%  to 
27%  over  the  ap  range  (Marcos,  et  ah,  2009).  Figure  17b  shows  CHAMP  vs  empirical 
models  (Marcos,  et  ah,  2009)  for  the  same  time  period,  but  with  no  hemisphere  or  local 
time  binning.  The  JB08  model  is  in  green;  NRLMSIS  and  JB06  are  purple  and  orange 
respectively.  JB08  and  NRLMSIS  means  are  close  to  unity  for  all  bins,  but  JB06  shows  a 
decrease  of  about  20%.  Standard  deviations  for  JB08  increase  from  about  16%  (lowest  ap 
bin)  to  about  27%  (highest  ap  bin)  and  are  therefore  about  9%  lower  for  the  lowest  bin 
and  only  four  percent  lower  for  high  geomagnetic  activity.  Both  JB08  and  NRLMSIS 
have  standard  deviations  of  about  32%  in  the  highest  bin.  It  is  noted  again  that  each  ap 
bin  in  this  figure  contains  data  for  all  latitudes  and  local  times;  finer  binning  could  reduce 
errors. 

Figure  18a  shows  CHAMP  data  in  10  degree  latitude  bins  (data  for  each  bin  plotted  at 
start  of  interval)  for  day,  left  side  and  night,  right  side.  Mean  values  are  at  the  top  and 
standard  deviations  are  at  the  bottom.  For  each  chart  data  are  plotted  for  the  lowest  ap 
bin  (blue)  and  ap>100  (red).  Means  are  fairly  constant  vs  latitude  on  the  dayside  and 
increase  by  about  15%  for  high  ap  on  the  nightside.  Nightside  errors  are  higher  than  those 
on  the  dayside  for  both  ap  bins.  Day  errors  increase  from  about  24%  to  36%  whereas  the 
night  values  increase  from  about  29%  to  37%.  Figure  18b  shows  CHAMP  vs  JB08  means 
(top)  and  standard  deviations  (bottom).  Mean  values  tend  to  increase  by  about  10%  from 
low  to  high  ap.  Standard  deviations  increase  from  about  15%  to  19%  at  low  ap  with 
increasing  latitude.  For  high  ap,  errors  increase  from  about  25%  to  37%  for  ap  100-200 
but  drop  to  19%  for  the  highest  bin,  possibly  reflecting  impact  of  the  new  geomagnetic 
storm  algorithm  implemented  in  JB08  (Burke,  et  ah,  2007).  Note  that  at  high  southern 
hemispheres  the  error  is  -26%.  Thus  TIEGCM  errors  in  the  100-200  ap  bin  are  similar  to 
those  for  JB08  the  northern  hemisphere.  For  the  ap  >200  case,  JB08  errors  are  lower. 

SETA-2  to  TIEGCM  data  in  Figure  19a  also  show  mean  ratios  at  top  and  standard 
deviation  at  bottom.  Blue  points  are  daytime  data  and  red  are  for  nighttime  bins.  Mean 
ratios  in  the  northern  hemisphere  increase  on  average  from  about  0.96  to  1.2  as  ap 
increases  from  lowest  to  highest  ap  bin  and  standard  deviations  rise  from  1 1  to  18%. 
These  lower  standard  deviations  are  assumed  to  be  mainly  an  altitude  effect,  consistent 
with  one-day  average  drag  data  showing  that  the  standard  deviation  increases  with 
altitude  from  about  10%  at  220  km  to  -15%  at  400  km  (Marcos,  et  ah,  2006).  Local  time 


Approved  for  public  release;  distribution  is  unlimited. 


16 


variations  are  small;  daytime  errors  are  a  few  percent  higher  in  the  two  highest  ap  bins. 
Corresponding  JB08  (green  data)  means  for  day  local  times  in  the  northern  hemisphere  in 
Figure  19  were  fairly  constant  varying  from  about  0.94  to  1.02  over  the  ap  span  (Lin,  et  ah, 
2011).  Standard  deviations  varied  from  16  to  21%  from  lowest  to  highest  bin  values. 
Increases  in  error  with  increases  in  ap  bin  for  JB08  (green  points)  are  very  similar  to 
those  of  J70  (red  squares).  There  is  not  an  improvement  for  JB08  in  the  higher  ap  bins  as 
was  observed  for  CHAMP. 

Figure  20a  shows  SETA-2  data  in  10  degree  latitude  bins  (data  for  each  bin  plotted  at 
start  of  interval)  for  day,  left  side,  and  night,  right  side.  Mean  values  are  at  the  top  and 
standard  deviations  are  at  the  bottom.  For  each  chart  data  are  plotted  for  the  lowest  ap 
bin  (blue)  and  ap>100  (red).  The  night  values  are  higher  than  those  for  day  in  all  cases. 

The  difference  may  be  due  to  modeling  or  to  a  bias  in  the  SETA  processing.  Reference  to 
Figure  5  shows  an  average  dayside  altitude  of  ~  185  km  and  an  average  nightside  altitude 
-230  km.  In  the  high  ap  bin,  standard  deviations  from  equator  to  pole  increase  for 
TIEGCM  ratios  from  14  to  28%  (day)  and  19  to  26%  (night).  Ratios  to  JB08  are  in  Figure 
20b  for  three  models.  The  JB08  errors  for  the  highest  ap  bin  (left  figures)  increase 
froml  1%  to  20%.  Thus  the  JB08  model  errors  vs  latitude  are  smaller  overall. 

8.  Summary 

General  Circulation  Models  have  been  a  critical  research  tool  for  interpreting  data  and 
advancing  understanding  physics  of  the  thermosphere  and  ionosphere.  Simulations 
depend  on  the  model  inputs  and  boundary  conditions,  including,  as  noted  in  Section  lc, 
solar  radiation,  magnetospheric  currents,  auroral  particle  fluxes,  and  global-scale 
atmospheric  waves  from  below.  There  are  large  uncertainties  in  these  highly  variable 
inputs  at  any  given  time.  Observations  are  relatively  sparse,  and  temporaFspatial 
extrapolation  of  the  observations  to  adequately  specify  boundary  conditions  is  difficult. 
Validation  efforts  detennine  fidelity  of  models  and  guide  further  development.  This 
preliminary  assessment  of  TIEGCM  version  1.9.1  using  global  neutral  density 
measurements  from  satellite  accelerometer  data  has  shown: 

a.  Solar  cycle  variations:  TIEGCM  overestimates  density  during  periods  near  solar 
activity  minimum,  as  do  empirical  models  since  in  both  cases  the  solar  input  was 
parameterized  by  the  F10  index.  The  model  has  the  capability  to  input  high  resolution 
solar  EUV  data.  Implementation  of  selected  lines  has  shown  promise  to  improve 
specification  but  the  wavelengths  and  their  relative  importance  need  further  study. 

Detailed  forthcoming  solar  irradiance  measurements  from  the  EVE  experiment  on  SDO 
will  eventually  permit  the  full  capability  of  physical  models  to  represent  thermosphere 
response  to  variations  in  solar  EUV.  We  note  that  EUV  bursts  from  solar  flares,  not 
described  in  this  Memo,  are  another  source  of  abrupt  density  perturbations  (Sutton,  et  al., 
2006).  Qian  and  Solomon,  2012  have  shown  that  TIEGM  has  the  capability  to  handle 
these  given  the  correct  inputs. 

b.  Semiannual  variations:  The  physics  for  this  effect  did  not  fall  out  of  V  1.9.1  and  earlier 
versions  of  this  model.  An  empirical  seasonally  dependent  eddy  diffusion  coefficient  at 
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90  km  deduced  from  satellite  drag  data  near  400  km  was  invoked  to  give  a  realistic 
TIEGCM  response  in  version  1.9.4.  A  theoretical  solution  using  gravity  waves  and  tides 
is  being  sought.  The  trend  of  long  term  semiannual  variations  can  be  modeled  in  an 
operational  forecast  scheme.  If  irregular  contributions  that  perturb  this  trend  can  be 
accounted  for,  then  physical  models  would  exceed  the  capability  of  empirical  models  for 
density  semiannual  forecasts. 

c.  Local  time  variations:  Nightside  (18-06  hrs)  TIEGCM  model  to  data  ratios  are  on 
average  -15%  higher  than  on  the  dayside.  Data  assimilation  suggested  that  the  physical 
model  cooling  rates  need  adjustment.  The  feature  was  also  found  in  TIEGCM  model  V 
1 .9.4.  This  effect  requires  further  study. 

d.  Geomagnetic  storms:  TIEGCM  simulations  tended  to  underestimate  storm  response  at 
low  and  high  thermospheric  altitudes  due  to  lack  of  accuracy  of  heating  drivers,  mainly 
the  quality  of  input  electric  fields.  On-going  data  assimilation  efforts  show  that  can 
geomagnetic  storm  density  errors  (as  well  as  those  due  to  other  effects)  significantly 
reduce  these  errors. 

e.  The  current  statistical  accuracy  of  TIEGCM  is  approaching  that  of  empirical  models 
from  CHAMP  altitudes  down  to  the  reentry  region.  This  very  encouraging  finding  is 
consistent  with  validation  efforts  of  other  physical  models  at  CHAMP  altitudes  outside 
the  scope  of  this  memo  (e.g.  Codrescu,  et  ah,  2012;  Shim,  et  ah,  2011). 

The  ultimate  goal  of  satellite  drag  research  is  to  implement  an  accurate  assimilative, 
operational,  physics  based  forecast  model  with  5%  accuracy  and  3-5  day  forecasts.  With 
much  new  understanding  of  the  atmosphere  being  developed,  physics-based  modeling 
tools  currently  are  the  future  tool  for  density  specification  and  prediction.  Physics-based 
models  offer  an  advantage  by  providing  more  detailed,  time-dependent  structure  not 
available  in  the  semi-empirical  models.  The  use  of  input  proxies  for  solar  UV  and  EUV 
radiation,  of  statistical  patterns  of  particle  precipitation  and  plasma  convection  for  the 
calculation  of  Joule  heating,  and  of  seasonal  averages  for  the  amplitude  and  phase  of 
waves  propagating  from  below  leads  to  large  uncertainties  in  model  predictions. 
Extensive  research  on  the  required  fundamental  heating  processes  resulting  from  the 
EUV  and  solar  wind  interaction  with  the  coupled  magnetosphere-ionosphere- 
thermosphere  system  continues.  Geomagnetic  storms  remain  the  biggest 
challenge.  Much  of  our  imperfect  understanding  is  due  to  lack  of  observational  data  on 
high  latitude  heating  mechanisms.  Joule  heating  is  the  dominant  heating  mechanism 
affecting  the  neutral  density  during  disturbed  conditions.  Knowing  its  magnitude  and 
spatial  distribution  is  important  for  the  accuracy  of  neutral  density  specification  and 
atmospheric  drag  detennination. 

The  challenges  facing  operational  physical  models  are  scientific,  requiring  new  space 
environment  theories  of  solar  phenomena,  their  propagation  through  the  interplanetary 
environment  and  their  interaction  with  the  Earth’s  thermosphere-ionosphere- 
magnetosphere  system  and  new  space  environment  measurements.  The  problems  for 
space  weather  predictions  are  analogous  to  those  addressed  by  tropospheric  weather  for 
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which  the  basic  physics  has  been  known  for  decades.  In  spite  of  knowing  how  the 
equations  work,  they  require  input,  hence  the  need  for  operational  weather 
satellites.  Observations  are  critical  to  unravel  the  basic  physical  mechanisms  to  build 
better  models  with  improved  predictions,  but  there  will  always  be  a  need  for  input  to  the 
models.  The  alternative  to  measuring  the  inputs  on  the  necessary  temporal  and  spatial 
scales  is  to  leverage  all  available  measurements  of  state  parameters  and  steer  the  model 
results  in  the  right  direction  using  a  data  assimilation  scheme. 

Results  from  the  AFOSR  MURI  NADIR  effort  headed  by  the  University  of  Colorado 
(see  University  of  Colorado  MURI  website  for  complete  description)  are  an  obvious 
resource  to  be  assimilated  into  operational  physical  modeling.  The  MURI  team  has 
developed  research  products  that  can  be  transitioned  into  operational  models.  Some 
transition  highlights  were  reviewed  by  Fuller-Rowell  and  Forbes,  2010.  The 
fundamental  next  step  after  the  MURI  is  to  incorporate  the  benefits  of  a  physical  model 
within  the  assimilation  process.  In  the  future,  solar  wind  and  magnetosphere  models  may 
be  able  to  predict  important  external  drivers  for  thermosphere-ionosphere  models,  so  a 
few  hours  to  a  few  days  forecast  of  the  Joule  heating  index  and  therefore  neutral  density 
would  be  feasible.  The  MURI  program  attempts  to  address  these  model  limitations 
through  scientific  analyses  of  solar  phenomena  and  available  satellite  drag 
data.  Initiating  the  physical  model  with  both  direct  in-situ  and  indirect  observations  in  a 
statistically  rigorous  manner  provides  a  practical  approach  for  representing  the  time- 
dependent  conditions  of  the  thermosphere.  Implementation  of  a  Kalman  filtering  data 
assimilation  technique  by  the  space  physics  community  has  proven  to  be  non-trivial.  It 
requires  detailed  understanding  of  the  physics  based  model  as  well  as  estimation  theory 
and  astrodynamics.  Thus,  this  development  requires  cooperative  efforts  of  the  space 
physics  and  astrodynamics  communities.  AFRL  has  demonstrated  application  of  the 
technique  reduces  errors  during  a  major  geomagnetic  storm. 

In  summary,  uncertainties  in  neutral  density  variations  have  been  the  major  error  source 
for  LEO  determination.  The  problem  is  being  vigorously  and  fruitfully  attacked  by 
numerous  space  weather  studies  including  physical  modeling  research,  data  assimilation 
schemes,  predictive  solar  and  geomagnetic  indices,  upward  propagating  waves  and  in-situ 
measurements.  While  there  is  still  a  lot  of  research  needed,  the  tools  are  finally 
becoming  available.  The  culmination  of  these  efforts  will  be  steady  progress  in  meeting 
the  evolving,  previously  unattainable,  stringent  requirements  for  operations  in  the  satellite 
drag  environment. 
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Appendix:  Figures 


Figure  1:  Thermosphere  Heating  and  Dynamics  Processes. 
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Figure  2:  Solar  electromagnetic  and  solar  wind  heating  into  thermosphere. 
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Figure  3:  Assimilative  Physics-Based  Neutral  Density  Models. 


Figure  4:  Solar  activity  and  spacecraft  orbital  prediction. 
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Figure  5:  Typical  orbital  parameters  for  the  SETA  missions. 
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Figure  6:  TIGCM  underestimation  of  density  at  solstices  for  (top) 
satellite  drag  measured  by  two  different  satellites  showing  model/data 
ratio  vs  time  and  (Bottom)  model  to  density  from  CHAMP 
accelerometer  in  2002. 
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Modified  T1E-GCM  Eddy  Diffusivity  Scheme 


TIE-GCM  ( 2002_plvmeanadj )  /  CHAMP  Comparison:  DAYSIDE 
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Figure  7:  Top:  Variable  eddy  diffusion  coefficient  imposed  on  TIEGCM.  Middle: 
Revised  model  densities  vs  CHAMP  density  data.  Bottom:  Revised  model/CHAMP 
ratios. 
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Figure  8a:  Solar  cycle  variation  of  TIEGCM  and  NRLMSIS  annual  mean  values. 
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Figure  8b:  Solar  cycle  variation  of  TIEGCM  and  NRLMSIS  standard  deviations. 
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Figure  8c:  Top:  Solar  spectral  radiance;  Bottom:  Ratio  of  solar  maximum  to  solar 
minimum  radiance  vs  wavelength. 
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Figure  9:  CHAMP  &  GRACE  perigee  altitude  (top)  and  solar  flux  (bottom)  vs  time. 
Blue  dashed  lines  indicate  periods  studied. 
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Figure  10a:  CHAMP  2001  geomagnetic  storm  responses  at  high  solar  flux.  Left  top 
and  middle  panels  are  for  period  starting  day  304;  right  top  and  middle  panels  are 
for  periods  starting  day  357.  Bottom  panels  are  respective  ap  values. 
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Figure  10b:  TIEGCM  simulations  of  two  storms  in  Figure  10a.  Storm  1  is  in  top  two 
panels;  storm  2  is  in  bottom  two  panels. 
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Figure  10c:  Ratio  of  CHAMP  to  TIEGCM  densities  for  2001  high  solar  flux  storms 
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Figure  11a:  CHAMP  (top  panels)  and  GRACE  (middle  panels)  response  to  same 
high  solar  flux  storm.  Bottom  panels  are  ap  values. 
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Figure  lib:  Ratio  of  CHAMP  and  GRACE  to  TIEGCM  densities  for  2002  high  solar 
flux  storm. 
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Figure  12a:  CHAMP  and  GRACE  2007  geomagnetic  storm  responses  at  low  solar 
flux.  Top:  CHAMP  data;  middle:  GRACE  data.  Bottom  panels  are  ap  values. 
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Figure  12b:  TIEGCM  simulations  of  2007  storm  at  CHAMP  (top)  and  GRACE 
(bottom)  altitudes. 
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Figure  12c:  Ratio  of  CHAMP  and  GRACE  to  TIEGCM  densities  for  2007  low  solar 
flux  storm. 
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Figure  13a:  CHAMP  and  GRACE  2008  geomagnetic  storm  responses  at  low  solar 
flux.  Top:  CHAMP  data;  middle:  GRACE  data.  Bottom  panels  are  ap  values. 
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Figure  13b:  Ratio  of  CHAMP  (top)  and  GRACE  (bottom)  to  TIEGCM  densities  for 
2008  low  solar  flux  storm 
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Figure  14a:  SETA-2  neutral  mass  density  at  200  km  in  latitude-time  coordinates  for 
3-8  Aug  1982  for  local  times  near  1030  hours  (top)  and  2230  hours.  Geomagnetic 
indices  ap  and  Dst  are  in  bottom  panel. 
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Figure  14b:  TIEGCM  simulation  of  neutral  densities  normalized  to  200  km  at 
SETA-2  locations  for  d.  216-220  (4-8  Aug  1982). 

Top:  Dayside  (1030  LT);  Middle:  Nightside  (2230  LT).  Bottom:  Solar  wind  inputs. 
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Figure  14c:  Ratio  of  SETA-2  to  TIEGCM  density  for  Aug  1982  storm  for  Local 
times  near  1030  hours  (left)  and  2230  hours  (right). 
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Northern  Hemisphere  Southern  Hemisphere 


Figure  15.  Orbit-averaged  density  and  models  vs  time  in  Aug  1982.  Black:  SETA-2 
measurements,  Red:  J70  model,  Green:  JB08  model.  Orange:  NRLMSIS  model  and 
Blue:  TIEGCM  model. 
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Figure  16a:  Density  during  part  of  day  294,  2004,  vs  unconstrained  models.  Top: 
mass  density  from  CHAMP  (red),  NRLMSIS  (black)  and  TIEGCM  results  using 
SEE  EUV(pink),  EUVAC  (blue  and  green  with  different  time  resolutions  of 
geomagnetic  activity  data  -  see  text).  Middle:  Satellite  latitude.  Bottom:  Satellite 
longitude  (black)  and  local  time  (blue). 
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Figure  16b:  Same  as  for  Figure  16a  but  with  TIEGCM  runs  constrained  by 
adjusting  rate  factors  to  better  match  the  density  data.  In  the  top  panel  the  blue 
curve  is  the  TIEGCM/SEE  run  with  the  NO  cooling  factor  adjusted;  the  green  curve 
adds  a  joule  heating  adjustment.  The  pink  and  black  curves  are  as  in  Figure  16a. 


Approved  for  public  release;  distribution  is  unlimited. 


45 


CHAMP 


DOY 


Plot  scale=  2e-12 
Data 

Model:  mschampdn.met 
Model:  sochampdn.met 
Model:  splpchampdn.met 
Model:  sp2pchampdn.met 


Figure  16c:  Data  assimilation  using  variable  NO  cooling  rate  and  joule  heating 
parameters  for  TIEGCM  Nov  2004  CHAMP  data  simulation.  Top  panel:  Mean 
daily  density  for  CHAMP  (red);  NRLMSIS  (black)  and  TIEGCM  runs  (i) 
unconstrained  (pink),  (ii)  with  adjusted  NO  cooling  (blue)  and  (iii)  adjusted  NO  and 
joule  heating  (green).  Second  panel  is  the  prediction  efficiency:  the  observational 
variance  minus  the  mean  square  model  error,  normalized  by  the  observational 
variance;  Third  panel  is  the  mean  model  error.  Fourth  panel  is  the  root  mean 
square  model  error. 
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Figure  16d:  Latitude-UT  difference  in  neutral  winds  at  240  km  between  two  model 
runs  (constrained  and  unconstrained)  during  an  interval  of  geomagnetic  activity 
(day  315,  2004).  The  top  panel  is  the  meridional  wind,  the  bottom  panel  the  zonal 
wind.  The  constrained  model  had  the  NO  cooling  factor  fitted  to  have  its  densities 
better  match  the  CHAMP  densities. 
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Figure  17a:  Mean  (top)  and  standard  deviations  (bottom)  of  CHAMP/TIEGCM 
northern  hemisphere  data,  in  six  ap  bins,  for  day  (blue)  and  night  (red)  local  time 
bins. 
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Figure  17b:  Mean  (top)  and  standard  deviations  (bottom)  of  CHAMP/Empirical 
Model  data  (no  local  time  or  hemisphere  binning) 
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Latitude 


Figure  18a:  CHAMP  to  TIEGCM  ratios  vs  latitude  for  northern  hemisphere  data  in 
two  ap  bins,  0-10  (blue)  and  >100  (red).  Left:  daytime  local  times;  Right  nighttime 
local  times 
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Figure  18b:  CHAMP  to  JB08  ratios  in  four  ap  bins.  Data  cover  all  latitudes  and  are 
for  all  local  times.  Mean  values  (top)  and  standard  deviations  (bottom). 
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Figure  19a:  Mean  (top)  and  standard  deviations  (bottom)  of  SETA-2/TIEGCM 
northern  hemisphere  data,  in  six  ap  bins,  for  day  (blue)  and  night  (red)  local  time 
bins. 
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Figure  19b:  Mean  (top)and  standard  deviations  (bottom)  of  SETA-2  to  empirical 
models,  day  local  times  and  northern  hemisphere  data,  in  six  ap.  Ratios  of  data  are 
to:  Blue:  NRLMSIS,  Red:  J70,  Green:  JB08. 
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Figure  20a:  SETA-2  to  TIEGCM  ratios  vs  latitude  for  northern  hemisphere  data  in 
two  ap  bins,  0-10  (blue)  and  >100  (red).  Left:  daytime  local  times;  Right:  nighttime 
local  times 
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Figure  20b:  Northern  hemisphere  dayside  SETA-2  mean  values  (top)  and  standard 
deviations  (bottom)  for  data-to-empirical  model  ratios  binned  by  dayside  (left)  and 
nightside  (right).  Data  are  for  ap  0-10  (left)  and  ap  100-400  (right).  Ratios  are  of 
data  to:  Blue:  NRLMSIS,  Red:  J70,  Green:  JB08. 
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